Namibian hake model update, 2025

Author

John Kathena, Jim Ianelli

Published

August 16, 2025

Show the code
library(knitr)
hook_output <- knit_hooks$get("output")
knit_hooks$set(output = function(x, options) {
  # this hook is used only when the linewidth option is not NULL
  if (!is.null(n <- options$linewidth)) {
    x <- knitr:::split_lines(x)
    # any lines wider than n should be wrapped
    if (any(nchar(x) > n)) x <- strwrap(x, width = n)
    x <- paste(x, collapse = "\n")
  }
  hook_output(x, options)
})

knitr::opts_chunk$set(collapse = TRUE, comment = "  ", fig.align = "center", cache = FALSE, tidy.opts = list(width.cutoff = 80), tidy = TRUE)
knitr::opts_knit$set(root.dir = here::here())
# knitr::opts_chunk$set(warning=F, message=F, echo=F, results=F,fig.width=6, fig.height=5)

Assessment model runs

The original base-case model was evaluated for a number of features and extensions. These included focus on what data components were fit well and how improvements in consistency can be made. For the latter part, we found that the fits to the index and CPUE data were particularly poor and could be improved. This document is an update from the base-case model presented in 2024.

Model descriptions

The following table was developed based on testing the model with different assumptions and data sources. Key differences from the 2023 assessment configuration was the assumption that model estimation of variance terms was appropriate. This feature resulted in unacceptable residual patterns and essentially a complete down weighting of the index data. We used the assumed variance terms (CVs) for the indices in all of the following model configurations:

Model Description
Pre 2024 model configuration As specified in past assessments, estimated steepness and all variance terms
Base case (m0) Model with survey “minus group” to be ages 0, and 1 instead of 0, 1, and 2 as done in the past, steepness fixed at 0.7, q estimated, and time-varying fishery asymptotic selectivity specified.
m25 As base case but updated through 2025 .
m1 As 2024 base case but with survey catchability fixed at 1.0
m2 As 2024 base case but with survey catchability fixed at 0.5
m3 As 2024 base case but with natural mortality estimated
m4 As 2024 base case but with fishery selectivity allowed to be dome-shaped
m5 As 2024 base case but with stock-recruit steepness fixed at 0.5
m6 As 2024 base case but with stock-recruit steepness fixed at 0.9
Show the code
library(NamibianHake)
library(flextable)
library(here)
library(tidyverse)
library(ggridges)
theme_set(ggthemes::theme_few())
library(xtable)
library(kableExtra)
set_flextable_defaults(digits = 3, decimal.mark = ".", big.mark = ",", na_str = "<na>")
Show the code
mod_ref <- c("old_bc", "m25", "m0", "m1", "m2", "m3", "m4", "m5", "m6")
mod_dir <- c("old_bc", "m25", "m0", "m1", "m2", "m3", "m4", "m5", "m6")
mod_label <- c("2023 base case", "2025 Update", "2024 base case", "Model 1", "Model 2",
    "Model 3", "Model 4", "Model 5", "Model 6")


#---Main code that extracts all the results from the model lists above------
res <- get_results(mod_names. = mod_label, moddir = mod_dir)

modlst <- res$modlst
old_bc <- modlst[[1]]
m0 <- modlst[[3]]
m25 <- modlst[[2]]
moddiag <- res$moddiag
dfsrr <- data.frame()
for (i in 1:length(mod_ref)) {
    dfsrr <- rbind(dfsrr, data.frame(Model = mod_label[i], SSB = modlst[[i]]$SSB,
        R = modlst[[i]]$Pred_Rec))
}
mods <- data.frame()
for (i in 1:length(mod_ref)) {
    mods <- rbind(mods, data.frame(moddiag[[i]], Model = names(moddiag[i])))
}

Fits to index data

that the variance terms (CVs) were estimated. A more standard approach, where CVs are specified based on the data (e.g., from design-based sampling theory), was used and this performed better (Figure 1, Figure 2). In the previous assessment, the fit tot he early period of CPUE data was better, but in this case the CV was estimated to be about 10%, and extremely low value for this type of data (Figure 3).

Similarly, we re-evaluated the ability of this model to estimate the stock recruitment productivity parameter (steepness). It is extremely rare that sufficient data are available to freely estimate this, even with extensive high-quality index data. Therefore we evaluated what assumptions were taken elsewhere (for this and related species, e.g., South African hake) and ran the models with steepness fixed at 0.7.

Figure 1: Base case model fits to main survey data compared to the previous assessment.
Figure 2: Base case model fits to the CPUE index 3 data compared to the previous assessment.
Show the code
dfcpue <- rbind(data.frame(Year = 1964:2024, Obs = m0$Obs_CPUE_1, predicted = m0$e_CPUE_1,
    Model = "Base case 2024"), data.frame(Year = 1964:2025, Obs = m25$Obs_CPUE_1,
    predicted = m25$e_CPUE_1, Model = "Base case 2025"))
dfcpue |>
    filter(Obs > 0) |>
    ggplot(aes(x = Year, y = Obs, color = Model)) + geom_point(color = "black") +
    geom_line(aes(y = predicted)) + geom_point(aes(y = predicted)) + ggtitle("CPUE 1") +
    ylim(0, NA)
ggsave(here("mods", "figs", "cpue1.png"), width = 6, height = 5)
Figure 3: Base case model fits to the CPUE index 1 data compared to the previous assessment.
Show the code
dfcpue <- rbind(data.frame(Year = 1964:2024, Obs = m0$Obs_CPUE_3, predicted = m0$e_CPUE_3,
    Model = "Base"), data.frame(Year = 1964:2025, Obs = m25$Obs_CPUE_3, predicted = m25$e_CPUE_3,
    Model = "Base updated 2025"), data.frame(Year = 1964:2024, Obs = modlst[[3]]$Obs_CPUE_3,
    predicted = modlst[[3]]$e_CPUE_3, Model = "Model 1, q=1.0"), data.frame(Year = 1964:2024,
    Obs = modlst[[4]]$Obs_CPUE_3, predicted = modlst[[4]]$e_CPUE_3, Model = "Model 2, q=0.5"))
dfcpue |>
    filter(Obs > 0) |>
    ggplot(aes(x = Year, y = Obs, color = Model)) + geom_point(color = "black") +
    geom_line(aes(y = predicted)) + geom_point(aes(y = predicted)) + ggtitle("CPUE 1") +
    ylim(0, NA)
ggsave(here("mods", "figs", "cpue1_all.png"), width = 6, height = 5)
Figure 4: Base case model fits to the CPUE index 3 data compared to 2025 and also models 1 and 2.
Figure 5: Base case model fits to main survey data compared to models 1 and 2.

Age composition fits

The age composition fits for the base case model are shown inf Figure 6 and Figure 7. The base case model uses a ‘minus group’ equal to ‘1’ for the survey data and for the fishery it was set to 2 (as was done previously).

Show the code
PlotAgeFit(x = m25, title = "Base case 2025", type = "fishery", fage = 2, lage = 7) +
    ggthemes::theme_few(base_size = 10)
ggsave(here("mods", "figs", "Age_comp_fish.png"), width = 9, height = 8)
Figure 6: Base case model fits to the fishery age composition data. Note that the base case model uses a ‘minus group’ equal to ‘2’ for the fishery data.
Show the code
PlotAgeFit(x = m25, title = "Base case 2025", type = "survey1", fage = 1, lage = 7) +
    ggthemes::theme_few(base_size = 10)
ggsave(here("mods", "figs", "Age_comp_surv.png"), width = 9, height = 8)
Figure 7: Base case model fits to the survey age composition data. Note that the base case model uses a ‘minus group’ equal to ‘1’ for the survey data.

Selectivity

Recruitment estimates

The recruitment results are consistent with those seen for the spawning biomass. The recent models shows a decline in recruitment in general, and that given the new data, the estimates for the 2024 year-class dropped from the mean value assumed from the 2024 assessment (Figure 11). Compared to the base case model, the alternative 2024 models show a range of recruitment estimates (Figure 12).

Show the code
mods |>
    filter(Model %in% mod_label[3:2], Year > 2010, Variable == "R") |>
    ggplot(aes(x = Year, y = value, ymin = ymin, ymax = ymax, type = Model, fill = Model)) +
    geom_errorbar(width = 0.95, position = "dodge", alpha = 0.3) + geom_bar(width = 0.95,
    stat = "Identity", position = "dodge") + ggthemes::theme_few() + ylab("Recruitment age 0") +
    xlab("Year")
ggsave(here("mods", "figs", "rec.png"), width = 6, height = 5)
Figure 11: Alternative model results on recruitment estimates.
Show the code
mods |>
    filter(Model %in% mod_label[2:9], Year > 2010, Variable == "R") |>
    ggplot(aes(x = Year, y = value, ymin = ymin, ymax = ymax, type = Model, fill = Model)) +
    geom_errorbar(width = 0.95, position = "dodge", alpha = 0.3) + geom_bar(width = 0.95,
    stat = "Identity", position = "dodge") + ggthemes::theme_few() + ylab("Recruitment age 0") +
    xlab("Year")
ggsave(here("mods", "figs", "recalt.png"), width = 6, height = 5)
Figure 12: Alternative model results on recruitment estimates.

Stock recruitment relationships

For each of the models there were some specified and estimated differences in the stock-recruitment relationships (Figure 13). When overlain with the recruitment estimates, these relationships appear to be relatively poorly defined (Figure 14).

Show the code
p1 <- dfsrr |>
    ggplot(aes(x = SSB, y = R, color = Model, shape = Model)) + geom_point() + geom_line() +
    xlim(c(0, 4000))
p1
   Warning: The shape palette can deal with a maximum of 6 discrete values because more
   than 6 becomes difficult to discriminate
   ℹ you have requested 9 values. Consider specifying shapes manually if you need
     that many of them.
   Warning: Removed 67 rows containing missing values or values outside the scale range
   (`geom_point()`).
   Warning: Removed 7 rows containing missing values or values outside the scale range
   (`geom_line()`).
ggsave(here("mods", "figs", "srr_curves.png"), width = 6, height = 5)
   Warning: The shape palette can deal with a maximum of 6 discrete values because more
   than 6 becomes difficult to discriminate
   ℹ you have requested 9 values. Consider specifying shapes manually if you need
     that many of them.
   Warning: Removed 67 rows containing missing values or values outside the scale range
   (`geom_point()`).
   Warning: Removed 7 rows containing missing values or values outside the scale range
   (`geom_line()`).
Figure 13: Stock-recruitment curves noted for the previous base-case and the updated model alternatives.
Show the code
# |
mods |>
    filter(Model %in% mod_label[2:9], Variable == "R" | Variable == "SSB") |>
    select(c(1:3, 6)) |>
    pivot_wider(names_from = Variable, values_from = value) |>
    ggplot(aes(x = SSB, y = R, label = Year, color = Model, fill = Model)) + geom_text(alpha = 0.7,
    size = 3) + ggthemes::theme_few() + facet_wrap(. ~ Model) + geom_line(data = dfsrr |>
    filter(Model %in% mod_label[2:9]), aes(y = R, x = SSB, color = Model), inherit.aes = FALSE)
ggsave(here("mods", "figs", "srr_yrs.png"), width = 6, height = 5)
Figure 14: Stock-recruitment curves and year-class estimates for alternative models

To show the history relative to the replacement yield, we can plot the SSB/Bmsy and the catch/replacement yield (Figure 15). This shows a significant difference between the previous assessment and that proposed as the base-case for this year. The alternative models show a range of results, as shown in Figure 16.

Show the code
tmp <- mods |>
    select(Year, Model, Variable, value) |>
    filter(Model %in% mod_label[2:3], Year > 1980, Variable %in% c("Catch_RY", "B_Bmsy")) |>
    pivot_wider(names_from = Variable, values_from = value) |>
    arrange(Model, Year)  #|>
tmp |>
    ggplot(aes(x = B_Bmsy, label = Year, y = Catch_RY, shape = Model, color = Model,
        fill = Model)) + geom_path() + ggthemes::theme_few() + geom_point() + geom_text(alpha = 0.5) +
    geom_hline(yintercept = 1) + geom_vline(xintercept = 1) + xlab("SSB/Bmsy") +
    ylab("Catch / replacement Yield")
ggsave(here("mods", "figs", "kobe1.png"), width = 6, height = 5)
Figure 15: Base case model showing the relative SSB estimates compared to the 2024 assessment.
Show the code
tmp <- mods |>
    select(Year, Model, Variable, value) |>
    filter(Model %in% mod_label[2:3], Year > 1980, Variable %in% c("Catch_RY", "Depletion")) |>
    pivot_wider(names_from = Variable, values_from = value)
tmp <- mods |>
    select(Year, Model, Variable, value) |>
    filter(Model %in% mod_label[c(2:3, 7:8)], Year > 1980, Variable %in% c("Catch_RY",
        "B_Bmsy")) |>
    pivot_wider(names_from = Variable, values_from = value) |>
    arrange(Model, Year)  #|>
tmp |>
    ggplot(aes(x = B_Bmsy, label = Year, y = Catch_RY, shape = Model, color = Model,
        fill = Model)) + geom_path() + ggthemes::theme_few() + geom_point() + geom_text(alpha = 0.5) +
    geom_hline(yintercept = 1) + geom_vline(xintercept = 1) + xlab("SSB/Bmsy") +
    ylab("Catch / replacement Yield")
ggsave(here("mods", "figs", "kobe2.png"), width = 6, height = 5)
Figure 16: Base case model showing the relative SSB estimates compared to the model alternatives assessment.

Stock status comparisons

The Table 1 show the results of the base case model compared to the alternative models for a number of key statistics. Among these, the best fitting models were the base-case and Model 5. However, model 5 fits indicated that the basis for fitting the data was virtually identical and that the slight benefit based on the AIC arises from the stock-recruitment residuals. For advice and consistency with productivity assumptions (i.e., the specification of the steepness parameter fixed at 0.7 as done for South Africa) we selected the base-case scenario.

Show the code
# |
dftmp <- NULL
mod_scen <- c(2:3, 5:8)
filler <- " "
names(filler) <- "-ln(Likelihood)"
for (ii in mod_scen) {
    x <- modlst[[ii]]
    nll <- round(x$ObjFun, 0)
    names(nll) <- "Overall"
    CPUE <- round(x$CPUE_Like, 0)
    names(CPUE) <- "CPUE"
    surv <- round(x$Survey_Like, 0)
    names(surv) <- "Survey"
    caa <- round(x$CAA_Likelihood, 0)
    names(caa) <- "Commercial CAA"
    caas <- round(x$CAAS_Likelihood, 0)
    names(caas) <- "Survey CAA"
    oneyrold <- round(x$Oneyearold_Likelihood, 0)
    names(oneyrold) <- "One yr-old biomass"
    rec <- round(x$RecRes_Likelihood, 0)
    names(rec) <- "Rec. resids."
    datalike <- nll - rec
    names(datalike) <- "Data likelihood sub-total"
    NumPars <- round(x$Npars, 0)
    names(NumPars) <- "Number parameters"
    AIC <- round(x$Akaike, 1)
    names(AIC) <- "AIC"
    v <- c(filler, nll, CPUE, surv, caa, caas, oneyrold, datalike, rec, NumPars,
        AIC)
    dftmp <- cbind(dftmp, v)
}
dftmp <- data.frame(rownames(dftmp), dftmp, row.names = NULL)
names(dftmp) <- c("Component", mod_label[mod_scen])
tabcap <- "Model fits to data components"
flextable(dftmp) |>
    set_caption(caption = tabcap) |>
    colformat_double() |>
    autofit()
Table 1: Model fits to data components

Component

2025 Update

2024 base case

Model 2

Model 3

Model 4

Model 5

-ln(Likelihood)

Overall

-44

-42

30

-36

-44

-50

CPUE

-48

-47

-46

-45

-44

-47

Survey

-33

-32

-18

-31

-30

-32

Commercial CAA

-128

-124

-68

-126

-133

-124

Survey CAA

91

90

113

86

91

90

One yr-old biomass

68

68

73

69

67

68

Data likelihood sub-total

-60

-57

21

-52

-61

-57

Rec. resids.

16

15

9

16

17

7

Number parameters

80

79

79

78

89

79

AIC

71.7

73.1

218.3

83.5

89.4

57.4

Show the code
dftmp <- NULL
mod_scen <- c(2:3)
for (ii in mod_scen) {
    x <- modlst[[ii]]
    Ksp <- round(x$KspSTD, 0)
    names(Ksp) <- "Unfished spawning biomass"
    Kexp <- round(x$KexpSTD, 0)
    names(Kexp) <- "Unfished expl. biomass"
    steepness <- round(x$Steep, 3)
    names(steepness) <- "SRR steepness"
    Cur_B <- round(x$Bstd[length(x$Bstd)], 3)
    names(Cur_B) <- "2024 SSB"
    Bmsy <- round(x$Spmsy, 0)
    names(Bmsy) <- "SSB_msy"
    Cur_B0 <- round(x$Cur_B0, 3)
    names(Cur_B0) <- "Current SSB over unfished"
    Cur_Bmsy <- round(x$Cur_Bmsy, 3)
    names(Cur_Bmsy) <- paste0("Current SSB over Bmsy")
    MSY <- round(x$MSY, 0)
    names(MSY) <- "MSY"
    ry <- round(x$aveRY_last5, 0)
    names(ry) <- paste0("recent 5-yr average replacement yield")
    ry_cur <- round(x$aveRY_last5 * x$Cur_Bmsy, 3)
    names(ry_cur) <- "Recent RY x current SSB /Bmsy"
    v <- c(Ksp, Kexp, steepness, Cur_B, Bmsy, Cur_B0, Cur_Bmsy, MSY, ry, ry_cur)
    dftmp <- cbind(dftmp, v)
}
dftmp <- data.frame(rownames(dftmp), dftmp, row.names = NULL)
names(dftmp) <- c("Statistic", mod_label[mod_scen])
tabcap <- "Selected management measures from alternative models. "
flextable(dftmp) |>
    set_caption(caption = tabcap) |>
    colformat_double() |>
    autofit()
# align=paste0('lll',strrep('r',length(mod_scen+1)))) kable(tab,
# caption.placement = 'top', include.rownames = FALSE, sanitize.text.function =
# function(x){x}) print(tab)
Table 2: Namibian hake stock estimates by model alternative.

Statistic

2025 Update

2024 base case

Unfished spawning biomass

2,846.000

2,856.000

Unfished expl. biomass

2,479.000

2,488.000

SRR steepness

0.700

0.700

2024 SSB

490.867

551.008

SSB_msy

942.000

942.000

Current SSB over unfished

0.172

0.193

Current SSB over Bmsy

0.521

0.585

MSY

264.000

263.000

recent 5-yr average replacement yield

151.000

157.000

Recent RY x current SSB /Bmsy

78.577

91.588

A comparison among the models for stock status and

Figure 17: Some reference point comparisons among model runs. aveRY_90 is the average replacement yeild since 1990, aveRY is the average replacement yield over the most recent 5 years before the current year, Cur_90 is the current (terminal year) SSB over the estimate from 1990, Cur_B0 is over the unfished estimate, and Cur_Bmsy is the ratio of current SSB over the Bmsy estimate.

Control rule application

Modeling Namibian hake survey data by species

Fisheries stock assessments require data that are reliably collected and compiled. Secondarily, assessment models should be configured to match the assumptions associated with the observed data. To account for survey trends between the two species of hake we applied the estimated observation errors to a simple state-space random walk model. This approach has a number of options for how process-errors can be specified and estimated. The observation model applies the observation-error variances \((\sigma_{j,t}^2)\) for the \(j^{th}\) species in year \(t (x_{j,t} )\). The indices are fit to latent state variables, e.g., the underlying population trend \(ln(\hat{Z}_{j,t})\) as follows:

\(ln(Z_{j,t} ) = ln(\hat{Z}_{j,t} )+\epsilon_{j,t}\) where \(\epsilon_{j,t}∼N(0,\sigma_{j,t}^2 )\)

and the state equation and associated process error variance \(\sigma_{PE}^2\) is defined as

\(ln(\hat{Z}_{j,t+1} = ln(\hat{Z}_{j,t+} )+\eta_{j,t},\) where \(\eta_{j,t}∼N(0,\tau_j^2 ).\)

The process error variances \(\tau_j^2\) (which may or may not vary across indices) are fixed effect parameters and the unobserved species combined population \(ln(Z_{j,t} )\) is estimated as a series of random effects. The model is fit using maximum likelihood estimation in TMB using the R package “rema” (Sullivan 2022, rema library). The survey data for each species was used with CVs applied for observation error specifications. The values for \(\tau_j^2\) were tested for each species and found to be similar so they were set to the same values.

The above analysis provides a summary of the model runs and the design of a control rule that accounts for the signals in the data on the different species.

Application to the Namibian hake stocks

The control rule was first applied to the Namibian hake stocks using the survey data and the relative proportions of the two species. However, we wish to have a more robust control rule that accounts for the signals in the data on M. paradoxus biomass and have that be independent of M. capensis (e.g., ?@fig-relmean). For that case, we applied the survey smoothing model to paradoxus alone. The next step was to compute mean biomass over the period 1990-2024, and evaluate the adustment for different levels of \(\gamma\). The historical adjustments based on this aspect of the MP is shown in Figure 18.

Model runtime: 0.1 seconds stats::nlminb thinks the model has converged: mod\(opt\)convergence == 0 Maximum gradient component: 2.35e-09 Max gradient parameter: log_PE TMB:sdreport() was performed successfully for this model Historical biomass estimates for M. paradoxus and the mean biomass over the period 1990-2025.

Figure 18: Historical adjustments for different reactivities to M. paradoxus biomass.

Control rule developments

The situation for developing a two-species control rule where catches between species and the trend in overall biomass for both species is combined is challenging. For management purposes, the goal is to avoid incidental takes in the proportion of one species that exceeds the historical levels of depletion for either stock. Fortunately, survey data are available that can be used to distinguish trends in the relative biomass for both stocks. The design of the triggered control rule therefore must consider patterns in the relative biomass from the survey data, the absolute biomass of the combined catch and biomass as modeled from the combined-stocks assessment. This provides a pragmatic approach using available data.

The steps in the control rule would be first to run a simple model that projected the survey biomass and relative proportion of M. paradoxus. Then, given the mean proportion over the period, compute the adjustment needed to the overarching control rule for the management procedure. For example, the historical range based on the survey has been between about one third of the mean value (in the earliest part of the period) to about 70% above the mean proportion. This range (especially the lower value) was used as a semi-empirical way to develop a minimum stock size threshold as part of the control rule. That is, when the stock of M. paradoxus drops below 30% of the mean proportion of the combined stocks estimate, the TAC recommendation for the combined stock would be zero. So if proportion of M. paradoxus \((p_y)\) is greater than 30% of the long-term mean proportion, then

\(TAC_y=(\delta \sum_{y}^{y-4}\frac{RY_y} 5) \min(1.0,B_y/B_{MSY})^\lambda{0.5}^{-\lambda} )\min(1,\dot{p_y}/\dot{\bar{p}})^\gamma\)

where the rebuilding factor (\(\delta\)) is set to 0.8 when the spawning biomass is below \(B_{MSY}\) and 1.0 when above. In words, the TAC in year y is equal to the catch under the current control rule times the ratio of the spawning biomass relative to \(B_{MSY}\) (or proxy) and the externally estimated proportion of M. paradoxus from survey data. The second term relates \(B_{MSY}\) and is intended to take fast action (for \(\lambda >1.0\)) when the biomass falls below 0.5 of that value (a standard in many places to define “overfished”). The third term on the right hand side reflects the impact the M. paradoxus biomass projected from a survey smoother (described below) and adjusts the TAC advice downwards when the projected \(\dot{p_y}\) drops below the mean value. The values of \(\gamma, \lambda\) were evaluated and are shown in Table 3. We note that the specification of a survey linkage by the individual species provides an appropriate adjustment that reduces the exploitation rate and prevents potential for “the point of recruitment impairment” (PRI).

For the control rule as specified, the reactivity of the TAC advice to changes in either B/Bmsy or the biomass of M. paradoxus biomass can be adjusted. These are shown in Figure 19. For the purposes of this analysis, we set \(\gamma = .25\). For \(\lambda\), most models evaluated were above 0.5 of \(B_{MSY}\) so there was no added adjustment beyond the \(R=0.8\). So, following the TAC as specified, the base-case model was below \(B_{MSY}\) so \(\delta=0.8\) with the average replacement yield over the last 5 years is recommended, with the adjustment based on the relative biomass of M. paradoxus and the long-term mean biomass of the species.

The TAC advice for the combined stocks is then given in Table 3.

Show the code
# |
dftmp <- NULL
mod_scen <- c(2:3)
ii <- 2
names(cur_mp) <- "Ratio of paradoxus to mean "
for (ii in mod_scen) {
    x <- modlst[[ii]]
    Cur_Bmsy <- round(x$Cur_Bmsy, 3)
    names(Cur_Bmsy) <- paste0("Current B over Bmsy")
    ry <- round(x$aveRY_last5, 0)
    names(ry) <- paste0("recent 5-yr avg repl. yield")

    ry_cur <- round(x$aveRY_last5 * min(1, x$Cur_Bmsy), 3)
    # names(ry_cur) <- ('Recent RY x current B/Bmsy')

    lambda <- 1.5
    gamma <- 0.25
    if (Cur_Bmsy > 0.5)
        adj <- 1 else adj <- Cur_Bmsy^lambda/0.5^lambda
    tac1 <- round(0.8 * ry * adj * min(1, cur_mp)^gamma, 2)
    names(tac1) <- paste0("opt 1, lambda=", lambda, " gamma=", gamma)

    lambda <- 2
    gamma <- 0.25
    if (Cur_Bmsy > 0.5)
        adj <- 1 else adj <- Cur_Bmsy^lambda/0.5^lambda
    tac2 <- round(0.8 * ry * adj * min(1, cur_mp)^gamma, 2)
    names(tac2) <- paste0("opt 2, lambda=", lambda, " gamma=", gamma)

    lambda <- 2
    gamma <- 0.1
    if (Cur_Bmsy > 0.5)
        adj <- 1 else adj <- Cur_Bmsy^lambda/0.5^lambda
    adj
    Cur_Bmsy
    tac3 <- round(0.8 * ry * adj * min(1, cur_mp)^gamma, 2)
    names(tac3) <- paste0("opt 3, lambda=", lambda, " gamma=", gamma)

    lambda <- 2
    gamma <- 1
    if (Cur_Bmsy > 0.5)
        adj <- 1 else adj <- Cur_Bmsy^lambda/0.5^lambda
    adj
    Cur_Bmsy
    tac4 <- round(0.8 * ry * adj * min(1, cur_mp)^gamma, 2)
    names(tac4) <- paste0("opt 4, lambda=", lambda, " gamma=", gamma)

    filler <- ""
    names(filler) <- "TAC"

    v <- c(ry, Cur_Bmsy, cur_mp, filler, tac1, tac2, tac3, tac4)
    dftmp <- cbind(dftmp, v)
}
dftmp <- data.frame(rownames(dftmp), dftmp, row.names = NULL)
names(dftmp) <- c("Statistic", mod_label[mod_scen])
tabcap <- "Selected management measures from alternative models. "
flextable(dftmp) |>
    set_caption(caption = tabcap) |>
    colformat_double() |>
    hline(i = 3) |>
    autofit()
Table 3: Namibian hake specification of different options for TAC considerations

Statistic

2025 Update

2024 base case

recent 5-yr avg repl. yield

151

157

Current B over Bmsy

0.521

0.585

Ratio of paradoxus to mean

0.931

0.931

TAC

opt 1, lambda=1.5 gamma=0.25

118.66

123.37

opt 2, lambda=2 gamma=0.25

118.66

123.37

opt 3, lambda=2 gamma=0.1

119.94

124.71

opt 4, lambda=2 gamma=1

112.46

116.93

Show the code
#--- Plot the adjustment factors for different reactivities to M. paradoxus biomass----
df01 <- tibble(relbiom = 1:20/20, `Gamma = 0.1` = (1:20/20)^0.1, `Gamma = 0.5` = (1:20/20)^0.5,
    `Gamma = 0.25` = (1:20/20)^0.25, ) |>
    pivot_longer(cols = -relbiom, names_to = "Reactivity", values_to = "adjustment")
df01 |>
    ggplot(aes(y = adjustment, x = relbiom, color = Reactivity)) + geom_line() +
    xlab("Relative M. paradoxus biomass") + ylab("TAC adjustment") + ggthemes::theme_few()  #+ ylab('TAC adjustment based on M. paradoxus biomass')```{r lam_gam}
ggsave(here("mods", "figs", "gamma.png"), width = 6, height = 5)
Figure 19: TAC adjustments for different levels of current biomass of M. paradoxus biomass.
Show the code
#---Plot the adjustment factors for different reactivities to M. paradoxus biomass----
df01 <- tibble()
df01 <- rbind(df01, tibble(relbiom = 1:100/100, lambda = rep(0.5, 100)), tibble(relbiom = 1:100/100,
    lambda = rep(1, 100)), tibble(relbiom = 1:100/100, lambda = rep(1.5, 100)), tibble(relbiom = 1:100/100,
    lambda = rep(2, 100))) |>
    mutate(adjustment = if_else(relbiom < 0.5, relbiom^lambda/0.5^lambda, 1), lambda = as.factor(lambda))
df01 |>
    ggplot(aes(y = adjustment, x = relbiom, color = lambda)) + geom_line() + xlab("Spawning biomass relative Bmsy") +
    ylab("TAC adjustment") + ggthemes::theme_few()  #+ ylab('TAC adjustment based on M. paradoxus biomass')
ggsave(here("mods", "figs", "lambda.png"), width = 6, height = 5)
Figure 20: TAC adjustments for different levels of current biomass relative to Bmsy

Projections

A set of simple projections were established from the base case model. These were based on the 2024 catch estimates (lower and higher values) for contrast relative to depletion Figure 21.

Figure 21: 40-year projections at different levels of catch related to the catch in 2025 (e.g., 0.5*C_25 is constant catch at 50% of the catch estimated for 2025). The scenario labeled ‘Replacement yld’ is just the catch set to the replacement yield. The cases with ‘10yr’ represent catch scenarios where the catch is increased (and decreased) by 10% of the 2025 catch and thereafter set equal to the 2025 catch.